Εργαστηριακή Άσκηση 2

Σκοπός της δεύτερης σειράς ασκήσεων είναι η εξοικείωση με τις συναρτήσεις σχεδιασμού φίλτρων πεπερασμένης κρουστικής απόκρισης (FIR) και την υλοποίησή τους στο MATLAB. Προτού ξεκινήσετε την άσκηση θα πρέπει να μελετήσετε με προσοχή το Κεφάλαιο 1 και, ειδικότερα, την παράγραφο 1.3 του τεύχους του μαθήματος. Το MATLAB (www.mathworks.com) είναι ένα διαδραστικό εμπορικό πρόγραμμα (Windows, Linux, Unix) με το οποίο μπορείτε να κάνετε εύκολα αριθμητικές πράξεις με πίνακες. Μπορεί να το έχετε εγκατεστημένο τοπικά, στον προσωπικό σας υπολογιστή ή να εργάζεστε σε κάποιο Εργαστήριο Προσωπικών Υπολογιστών (ΕΠΥ) της Σχολής σας που διαθέτει το συγκεκριμένο λογισμικό.

Μέρος 1: Εισαγωγή

Στο MATLAB οι συναρτήσεις \(fft\) και \(ifft\) υποθέτουν ζεύγος μετασχηματισμού Fourier \(x(t)\) και \(X(f)\) υπολογισμένων σε μη αρνητικά διαστήματα \( t=[0:N-1]ts \) και \( f=[0:N-1]fo. \) Όπως έχετε ήδη δει στην Εργαστηριακή άσκηση 1, το άνω μισό μέρος του διαστήματος συχνοτήτων αντιστοιχεί στις αρνητικές συχνότητες του σήματος, όταν υπολογίζουμε το \(X(f)\) με τη βοήθεια της συνάρτησης \(fft\) (για σήματα πραγματικών τιμών, αυτά τα δυο μισά είναι κατοπτρικά ως προς το μέσο του διαστήματος). Ακριβώς το ίδιο ισχύει και για το άνω μισό μέρος του χρονικού διαστήματος, όταν το σήμα \(x(t)\) προκύπτει από τον αντίστροφο μετασχηματισμό Fourier μέσω της \(ifft\).

Το MATLAB διαθέτει τη συνάρτηση \(fftshift\) για να ολισθήσει κυκλικά τις τιμές του σήματος ή του μετασχηματισμού Fourier, ώστε να αντιστοιχούν σε κεντραρισμένα στο μηδέν αμφίπλευρα διαστήματα, δηλαδή, στις χρονικές στιγμές \(tb=[-ceil((N-1)/2): floor((N-1)/2)]ts\) ή στις συχνότητες \(fb=[–ceil((N-1)/2): floor((N-1)/2)]fo\). Με τον τρόπο αυτό μπορούμε να παράγουμε τα \(xb(t)\) και \(Xb(f)\) που αντιστοιχούν στην αμφίπλευρη αναπαράσταση του σήματος και του μετασχηματισμού Fourier.

Για να κατανοήσετε τα ανωτέρω θεωρείστε το διάνυσμα \([1 2 3 4]\) ως το αποτέλεσμα του \(FFT\) μήκους 4. Τότε, το πρώτο στοιχείο (1) είναι ο όρος dc, το τρίτο στοιχείο (3) είναι το σημείο στο μισό της συχνότητας δειγματοληψίας \(fs/2\), που μπορεί να εκληφθεί ότι αντιστοιχεί είτε στην \(–fs/2\) είτε στην \(+fs/2\). Τα στοιχεία 2 και 4 αντιστοιχούν στις συχνότητες \(+fs/4\) και \(–fs/4\). Εφαρμόζοντας την \(fftshift\), το στοιχείο 3 εμφανίζεται πρώτο, που σημαίνει ότι στο MATLAB αντιστοιχεί στην αρνητική συχνότητα \(–fs/2\), το επόμενο στοιχείο 4 αντιστοιχεί στη συχνότητα \(–fs/4\) ακολουθούμενο από το dc και τη συχνότητα \(+fs/4\). Για ένα μετασχηματισμό περιττού μήκους, δεν υφίσταται σημείο για το \(±fs/2\). Έτσι για το διάνυσμα \([1 2 3]\), η εφαρμογή της \(fftshift\) θα δώσει τα στοιχεία που αντιστοιχούν στις συχνότητες \(–fs/3, 0, +fs/3\).

Εκτός του ότι παράγουν εξόδους με τις αρνητικές συχνότητες ή χρόνους στο άνω μισό του διανύσματος, αμφότερες οι συναρτήσεις \(fft\) και \(ifft\) αναμένουν ως είσοδο διάνυσμα με την ίδια μορφή, αφού προφανώς ισχύουν οι ταυτότητες

\(h == \text{ifft}(\text{fft}(h))\) και \(H == \text{fft}(\text{ifft}(H))\)

Η πρώτη υποδεικνύει ότι η είσοδος της \(ifft\) πρέπει να είναι αντεστραμμένη, όπως την παράγει η \(fft\), και η δεύτερη ότι η είσοδος της \(fft\) πρέπει να είναι αντεστραμμένη, όπως την παράγει η \(ifft\).

Στο επόμενο σχήμα βλέπετε παραστατικά ένα ημιτονικό σήμα που έχει πολλαπλασιασθεί με παράθυρο Blackman τόσο στην αμφίπλευρη, όσο και την μονόπλευρη αναπαράστασή του.

import plotly.graph_objs as go
import numpy as np
from plotly.subplots import make_subplots


# Seed for reproducibility
np.random.seed(0)

# Generate a windowed random data
t = np.linspace(-40, 40, 80)  # 80 samples
window = np.zeros_like(t)
window[30:50] = 1  # non-zero between indices 30 and 50
random_data = np.random.randn(len(t))
x_t_windowed = window * random_data  # Apply window to random data

# Create a figure
fig = make_subplots()

# Add stem lines - one line plot for each stem
for i in range(len(t)):
    if window[i] != 0:  # Only plot stems for the windowed part
        fig.add_trace(go.Scatter(x=[t[i], t[i]], y=[0, x_t_windowed[i]], mode='lines', line=dict(color='blue'), showlegend=False))

# Add markers at the top of stems
fig.add_trace(go.Scatter(x=t, y=x_t_windowed, mode='markers', name='Windowed Random Data', marker=dict(color='blue')))

# Update yaxis and xaxis properties
fig.update_yaxes(title_text="Amplitude")
fig.update_xaxes(title_text="Time (samples)")

# Update layout
fig.update_layout(height=600, width=800, title_text="Windowed Random Data with Stem Lines")
fig.show(renderer='notebook')

import plotly.graph_objs as go
import numpy as np
from plotly.subplots import make_subplots


# Seed for reproducibility
np.random.seed(0)

# Generate windowed random data
t = np.linspace(-40, 40, 80)  # 80 samples
window = np.zeros_like(t)
window[30:50] = 1  # non-zero between indices 30 and 50
random_data = np.random.randn(len(t))
x_t_windowed = window * random_data  # Apply window to random data

# Μετατροπή των αρνητικών χρόνων στο τέλος
negative_time_indices = t < 0
positive_time_indices = t >= 0

# Συνδυασμός των δεδομένων με τα αρνητικά στο τέλος
t_reordered = np.concatenate((t[positive_time_indices], t[negative_time_indices]))
x_t_windowed_reordered = np.concatenate((x_t_windowed[positive_time_indices], x_t_windowed[negative_time_indices]))

# Create a figure
fig = make_subplots()

# Add stem lines for reordered data
for i in range(len(t_reordered)):
    fig.add_trace(go.Scatter(x=[t_reordered[i], t_reordered[i]], y=[0, x_t_windowed_reordered[i]], mode='lines', line=dict(color='blue'), showlegend=False))

# Add markers at the top of stems for reordered data
fig.add_trace(go.Scatter(x=t_reordered, y=x_t_windowed_reordered, mode='markers', name='Reordered Windowed Random Data', marker=dict(color='red')))

# Update yaxis and xaxis properties
fig.update_yaxes(title_text="Amplitude")
fig.update_xaxes(title_text="Time (samples)")

# Update layout
fig.update_layout(height=600, width=800, title_text="Reordered Windowed Random Data with Stem Lines")
fig.show(renderer='notebook')

Οι αντίστοιχοι μετασχηματισμοί Fourier, σε αμφίπλευρη και μονόπλευρη αναπαράσταση, φαίνονται στο επόμενο σχήμα:

Όταν τα \(x(t)\) και \(X(f)\) παράγονται από το MATLAB, δεν χρειάζεται κάποια ιδιαίτερη προσοχή, πλην της κυκλικής ολίσθησης σε περίπτωση που θέλουμε π.χ. να σχεδιάσουμε το αμφίπλευρο φάσμα ή σήμα. Όταν όμως ένα εκ των \(x(t)\) ή \(X(f)\) ορίζεται από τον χρήστη απαιτείται περισσότερη προσοχή, διότι, συνήθως χρησιμοποιούνται τα αμφίπλευρα σήματα ή φάσματα. Μπορείτε να μεταβείτε από τη μία αναπαράσταση στην άλλη ως εξής:

\(x = ifftshift(xb)\), \(X = fft(x)\), \(Xb = fftshift(X)\), εάν ξεκινάτε από αμφίπλευρο σήμα και θέλετε να καταλήξετε σε αμφίπλευρο φάσμα, και

\(X = ifftshift(Xb)\), \(x = ifft(X)\), \(xb = fftshift(x)\), εάν ξεκινάτε από αμφίπλευρο φάσμα και θέλετε να καταλήξετε σε αμφίπλευρο σήμα,

όπου η συνάρτηση \(ifftshift\) του MATLAB εκτελεί την αντίστροφη λειτουργία της \(fftshift\). Όταν το \(N\) είναι άρτιο, οι \(fftshift\) και \(ifftshift\) δίνουν το ίδιο αποτέλεσμα. Όταν όμως το \(N\) είναι περιττό αυτό δεν ισχύει και χρειάζεται προσοχή στη χρήση τους. Στην πράξη, η προσεκτική εφαρμογή των ανωτέρω έχει σημασία όταν υπολογίζεται η φάση του φάσματος. Το πλάτος του φάσματος δεν επηρεάζεται από την κυκλική ολίσθηση των στοιχείων που προκαλούν οι \(fftshift\) και \(ifftshift\) (δείτε ιδιότητες DFT).

Εξάσκηση

Δοκιμάστε στο παράθυρο εντολών τα ακόλουθα προκειμένου να εμπεδώσετε τη χρήση των συναρτήσεων \(fftshift\) και \(ifftshift\).

Show code cell source Hide code cell source

import numpy as np

# Create the array X
X = np.arange(-2, 3)  # equivalent to MATLAB's -2:2

# Apply fftshift and ifftshift
fftshifted_X = np.fft.fftshift(X)
ifftshifted_X = np.fft.ifftshift(X)

# Double fftshift and a combination of fftshift and ifftshift
Y = np.fft.fftshift(np.fft.fftshift(X))
Z = np.fft.ifftshift(np.fft.fftshift(X))

# Check if arrays are equal
X_equals_Y = np.array_equal(X, Y)
X_equals_Z = np.array_equal(X, Z)

print(f"X equals Y: {X_equals_Y}")
print(f"X equals Z: {X_equals_Z}")

X equals Y: False
X equals Z: True

Αντίστοιχος κώδικας MatLab

X = [-2:2]
fftshift(X)
ifftshift(X)
Y = fftshift(fftshift(X));
Z = ifftshift(fftshift(X));
isequal(X,Y)
isequal(X,Z)
        

Ερώτηση 1: Ποιο εκ των διανυσμάτων Υ και Ζ ισούται με το X; Γράψτε την απάντησή σας σε ένα αρχείο κειμένου lab2_nnnnn.txt, όπου nnnnn τα πέντε τελευταία νούμερα του αριθμού μητρώου σας, χρησιμοποιώντας το Notepad από το μενού των Windows (Start → Programs → Accessories → Notepad) και αποθηκεύστε το στον φάκελο My Documents. Θα υποβάλετε το αρχείο αυτό ηλεκτρονικά στο τέλος, αφού απαντήσετε και τις επόμενες ερωτήσεις, οπότε μπορείτε να τα αφήσετε ανοικτό. Ερώτηση 2: Επαναλάβατε με \(X=[-1:2]\). Τι παρατηρείτε; Γράψτε την απάντησή σας στο αρχείο κειμένου lab2_nnnnn.txt. Δοκιμάστε στο παράθυρο εντολών τα ακόλουθα δύο παραδείγματα για να εμπεδώσετε τη χρήση των συναρτήσεων \(fftshift\) και \(ifftshift\) σε συνδυασμό με τις \(fft\) και \(ifft\).

Show code cell source Hide code cell source

import numpy as np
from scipy.fft import fft, ifft, fftshift, ifftshift
import matplotlib.pyplot as plt
from ipywidgets import interact, Layout
import ipywidgets as widgets

# First part: Original signal and its FFT
xb = np.array([1, 2, 3, 4, 5, 4, 3, 2, 1])
x = ifftshift(xb)
X = fft(x)
Xb = fftshift(X)  # Spectrum with DC component in the center

# Second part: Low-pass filter effect
Xb_low_pass = np.array([0, 0, 1, 1, 1, 1, 1, 0, 0])
X_low_pass = ifftshift(Xb_low_pass)
x_low_pass = ifft(X_low_pass)
xb_low_pass = fftshift(x_low_pass.real)

# Function to plot
def plot_signals():
    fig1, axs = plt.subplots(2, 2, figsize=(15, 10))

    time_range = np.arange(-4, 5)
    
    axs[0, 0].stem(time_range, xb, linefmt='blue', markerfmt='bo', basefmt='r-')
    axs[0, 0].set_title('Original Signal')
    axs[0, 0].set_xlabel('Time (s)')
    axs[0, 0].set_ylabel('Amplitude')

    axs[0, 1].stem(time_range, np.abs(Xb), linefmt='red', markerfmt='ro', basefmt='r-')
    axs[0, 1].set_title('Magnitude Spectrum')
    axs[0, 1].set_xlabel('Frequency (Hz)')
    axs[0, 1].set_ylabel('Magnitude')

    axs[1, 0].stem(time_range, Xb_low_pass, linefmt='green', markerfmt='go', basefmt='r-')
    axs[1, 0].set_title('Low-pass Spectrum')
    axs[1, 0].set_xlabel('Frequency (Hz)')
    axs[1, 0].set_ylabel('Magnitude')

    axs[1, 1].stem(time_range, xb_low_pass, linefmt='purple', markerfmt='mo', basefmt='r-')
    axs[1, 1].set_title('Reconstructed Signal')
    axs[1, 1].set_xlabel('Time (s)')
    axs[1, 1].set_ylabel('Amplitude')

    plt.tight_layout()
    plt.show()

plot_signals()

Αντίστοιχος κώδικας MatLab

close all;
clear all;
clc;
xb=[1 2 3 4 5 4 3 2 1]                          % πραγματικό σήμα με άρτια συμμετρία
figure;
subplot (2,1,1);
plot([-4:4],xb);
ylabel('xb');
x=ifftshift(xb)                                 % το σήμα με τις αρνητικές συνιστώσες στο άνω μέρος
X=fft(x)                                        % FFT
Xb=fftshift(X)                                  % το φάσμα με τη dc συνιστώσα στο κέντρο, πραγματικές
                                                % τιμές με άρτια συμμετρία όπως αναμένεται
subplot (2,1,2);
plot([-4:4],Xb);
ylabel('Xb');
close all;
clear all;
clc;
Xb=[0 0 1 1 1 1 1 0 0]                          % φάσμα βαθυπερατού σήματος με άρτια συμμετρία
figure;
subplot (2,1,1);
plot([-4:4],Xb);
ylabel('Xb');
X=ifftshift(Xb)                                 % το φάσμα με τις αρνητικές συνιστώσες στο άνω μέρος
x=ifft(X)                                       % IFFT
xb=fftshift(x)                                  % πραγματικό σήμα με άρτια συμμετρία όπως αναμένεται
subplot (2,1,2);
plot([-4:4],xb);
ylabel('xb');
        

Ερώτηση 3: Τροποποιείστε το προηγούμενο παράδειγμα ώστε να ξεκινήσετε απευθείας με τον ορισμό του φάσματος του βαθυπερατού σήματος \(X\) όπως το αναμένει η \(ifft\). Γράψτε την απάντησή σας στο αρχείο κειμένου lab2_nnnnn.txt.

Show code cell source Hide code cell source

from scipy import signal
import scipy.io.wavfile
# Ανατρέξτε στην τεκμηρίωση της βιβλιοθήκης scipy.signal
# https://docs.scipy.org/doc/scipy/reference/signal.html
from scipy.fft import fft, fftfreq
from dash import Dash, dcc, html, Input, Output
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from numpy import random
import pandas
import plotly
import plotly.express as px
import plotly.graph_objs as go
import dash_bootstrap_components as dbc
import warnings
warnings.filterwarnings('ignore')
print("Libraries added successfully!")

Libraries added successfully!

Μέρος 2: Σχεδιασμός και υλοποίηση φίλτρων

Θα ασχοληθούμε με το Παράδειγμα 1.2 της παραγράφου 1.5 του τεύχους Μαθήματος. Το παράδειγμα αυτό παρουσιάζει δύο εναλλακτικές μεθόδους σχεδιασμού FIR φίλτρων: α) τη μέθοδο των παραθύρων και β) τη μέθοδο των ισοϋψών κυματώσεων τις οποίες εφαρμόζει στην περίπτωση βαθυπερατών φίλτρων.

Στο παράδειγμα, τα φίλτρα δοκιμάζονται σε ένα πραγματικό σήμα, s, το οποίο είναι αποθηκευμένο στο αρχείο sima.mat (binary αρχείο MATLAB). Πρόκειται για ένα σήμα sonar με φάσμα που εκτείνεται μέχρι περίπου τα 4 KHz και συχνότητα δειγματοληψίας Fs=8192 (είναι και αυτή αποθηκευμένη στο αρχείο sima.mat, μαζί με το σήμα).

Εδώ θα πειραματιστούμε με δύο σήματα: (i) το sonar του παραδείγματος, το οποίο εδώ διαβάζεται από ένα .txt αρχείο (έχει προέλθει με εξαγωγή του s από το MATLAB) και (ii) ένα σήμα μουσικής, το violin.wav (σήμα από μουσική βιολιού), το οποίο περιέχει υψηλότερες συχνότητες και έχει προέλθει με δειγματοληψία στα Fs_viol=44100 Hz.

Σήμα sonar

Show code cell source Hide code cell source

# Ανάγνωση δειγμάτων σήματος από txt file
with open('files/sima.txt') as f:
    s = [float(x) for x in f]
s=np.array(s)   
print('μέγεθος σήματος=', s.shape)
Fs=8192

μέγεθος σήματος= (6565,)

Στο πεδίο του χρόνου

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import FloatRangeSlider, VBox, HTML, Output, Layout
import ipywidgets as widgets
from IPython.display import display, clear_output

t = np.arange(0, len(s)) / Fs

# Plotting function
def plot_signal(t_range):
    plt.close('all')  # Close all existing figures
    fig = plt.figure(figsize=(10, 4))
    plt.plot(t, s, label='Signal', color='#00CC96')
    plt.xlim(t_range)
    plt.ylim([-0.05, 0.05])
    plt.title('Time Domain Plot of x')
    plt.xlabel('t (sec)')
    plt.ylabel('Amplitude')
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    return fig



# Create a range slider for time range selection
time_range_slider = FloatRangeSlider(
    value=[0, 0.8],
    min=0,
    max=t[-1],
    step=0.01,
    description='Time Range:',
    style={'description_width': 'initial'},
    layout=Layout(width='90%'),
    continuous_update=False
)

# Output widget for the plots
graph_output = widgets.Output()

# Function to update the plot inside the output widget based on slider changes
def response(change):
    with graph_output:
        clear_output(wait=True)  # Clear the previous plot before drawing a new one
        fig = plot_signal(time_range_slider.value)

# Observe changes in the slider and update the plot accordingly
time_range_slider.observe(response, names='value')

# Initial call to display the plot
response(None)


# Style the HTML element
html_label = widgets.HTML(
    value="""
        <h2 style='font-weight: bold; font-size: 30px; text-align: center;'>Signal and Fourier Transformation</h2>
    """
)

# Display the slider and the outputs
vbox_layout = Layout(display='flex', flex_flow='column', align_items='center', width='100%')
ui = widgets.VBox([html_label,
                   time_range_slider,
                   graph_output], layout=vbox_layout)

clear_output(wait=True)  # Clear the previous plot
display(ui)

Ακούμε το σήμα

Show code cell source Hide code cell source

import ipywidgets as widgets
from IPython.display import Audio, display

# Assuming `s` is your audio signal and `Fs` is the sampling rate
# If `s` is not defined, you need to define it or load an audio signal.

# Function to play the sound
def play_sound(b):
    display(Audio(data=20*s, rate=Fs))

# Create a button widget
button1 = widgets.Button(description="Play Sound")

# Display the button
display(button1)

# On button click, play the sound
button1.on_click(play_sound)

Φάσμα (spectrum)

Show code cell source Hide code cell source

import matplotlib.pyplot as plt
from scipy import signal

# Assuming `s` (your signal) and `Fs` (your sampling frequency) are defined
f, Pxx_den = signal.welch(s, Fs, noverlap=128, nperseg=256)

# Create matplotlib figure and axis
fig, ax = plt.subplots(figsize=(10, 6))

# Plot the Power Spectral Density (PSD) with a logarithmic scale on the y-axis
ax.semilogy(f, Pxx_den, color='#1F77B4')

# Update layout for a semilog plot similar to the Plotly configuration
ax.set_title('Φάσμα σήματος sonar')
ax.set_xlabel('Συχνότητα (Hz)')
ax.set_ylabel('Πυκνότητα φάσματος ισχύος')
ax.grid(True, which="both", linestyle='--', linewidth=0.5, color='gray')

# Display the figure
plt.tight_layout()  # Adjust layout to make room for the plot's elements
plt.show()

Σήμα βιολιού

Show code cell source Hide code cell source

f=open("files/violin.wav", 'rb')
Fs_viol, s_viol = scipy.io.wavfile.read(f)
print('Fs_viol=',Fs_viol, ' number of samples=',len(s_viol))
f.close()

Fs_viol= 44100  number of samples= 220500

Στο πεδίο του χρόνου

Show code cell source Hide code cell source

import matplotlib.pyplot as plt
import numpy as np

# Assuming s_viol (your signal) and Fs_viol (your sampling frequency) are defined
tvl = np.arange(0, len(s_viol)) / Fs_viol

# Create matplotlib figure and axis
fig, ax = plt.subplots(figsize=(10, 6))

# Plot the signal
ax.plot(tvl, s_viol, color='#1F77B4', linestyle='-')

# Update layout
ax.set_title('Signal over time')
ax.set_xlabel('Time (s)')
ax.set_ylabel('Amplitude')
ax.grid(True, which="both", linestyle='--', linewidth=0.5, color='gray')

# Display the figure
plt.tight_layout()  # Adjust layout to make room for the plot's elements
plt.show()

Show code cell source Hide code cell source

import ipywidgets as widgets
from IPython.display import display
import sounddevice as sd
import threading

# Define your sound and its sampling rate
# s_viol and Fs_viol should be defined before this point

def play_sound():
    sd.play(s_viol, Fs_viol)

def on_button_clicked(b):
    threading.Thread(target=play_sound).start()

# Create a button widget
play_button = widgets.Button(description="Play Sound")

# Link the button to the function that plays the sound
play_button.on_click(on_button_clicked)

# Display the button
display(play_button)

Φάσμα (spectrum) και Φασματόγραμμα (spectorgram)

Show code cell source Hide code cell source

import matplotlib.pyplot as plt
from scipy import signal
import numpy as np

# Assuming s_viol (your signal) and Fs_viol (your sampling frequency) are defined
f, Pxx_den = signal.welch(s_viol, Fs_viol, nperseg=1024, noverlap=256)

# Create matplotlib figure and axis
fig, ax = plt.subplots(figsize=(10, 6))

# Plot the Power Spectral Density (PSD) with a logarithmic scale on the y-axis
ax.semilogy(f, Pxx_den, color='#1F77B4')

# Update layout for a semilog plot
ax.set_title('Power Spectral Density')
ax.set_xlabel('Συχνότητα [Hz]')
ax.set_ylabel('Πυκνότητα φάσματος ισχύος [V^2/Hz]')
ax.grid(True, which="both", linestyle='--', linewidth=0.5, color='gray')


# Display the figure
plt.tight_layout()  # Adjust layout to make room for the plot's elements
plt.show()

Show code cell source Hide code cell source

import matplotlib.pyplot as plt
from scipy import signal
import numpy as np

# Assuming `s` (your signal) and `Fs` (your sampling frequency) are defined
f, tsp, Sxx = signal.spectrogram(s, Fs)

# Create matplotlib figure and axis
fig, ax = plt.subplots(figsize=(10, 6))

# Convert power to dB
Sxx_dB = 10 * np.log10(Sxx)

# Create the spectrogram
c = ax.pcolormesh(tsp, f, Sxx_dB)

# Add a colorbar
fig.colorbar(c, ax=ax, label='Intensity [dB]')

# Update layout
ax.set_title('Spectrogram')
ax.set_xlabel('Χρόνος [sec]')
ax.set_ylabel('Συχνότητα [Hz]')

# Display the figure
plt.show()

Βαθυπερατά φίλτρα

Η μέθοδος των παραθύρων

… …

Show code cell source Hide code cell source

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt

# Assuming Fs is defined somewhere in your code
Fs = 8000  # Example value, adjust as necessary
H = np.hstack((np.ones(int(Fs/8)), np.zeros(int(Fs-Fs/4)), np.ones(int(Fs/8))))

# Create the x values for the stem plot
x_values = np.arange(len(H))

# Create matplotlib figure and axis
fig, ax = plt.subplots(figsize=(10, 6))

# Create a stem plot
markerline, stemlines, baseline = ax.stem(x_values, H, linefmt='b-', markerfmt='bo', basefmt="r-")
plt.setp(stemlines, 'color', '#1F77B4')  # setting the color of the stem lines
plt.setp(markerline, 'color', '#1F77B4')  # setting the color of the markers

# Update layout
ax.set_title('Stem plot')
ax.set_xlabel('Index')
ax.set_ylabel('Amplitude')

# Show figure
plt.show()

Ορθογωνικό παράθυρο (απλή περικοπή της h)

Show code cell source Hide code cell source

%matplotlib inline
import ipywidgets as widgets
from IPython.display import display, clear_output
import matplotlib.pyplot as plt
import numpy as np


h = np.real(np.fft.ifft(H))
middle = int(len(h)/2)
h = np.hstack((h[middle:], h[:middle]))

h_variants = {
    'h32': h[middle-16:middle+16],
    'h64': h[middle-32:middle+32],
    'h128': h[middle-64:middle+64],
    'h160': h[middle-80:middle+80],
    'h256': h[middle-128:middle+128],
}

# Function to plot
def plot_stem(h_key):
    h_data = h_variants[h_key]
    x_values = np.arange(len(h_data))
    plt.close('all')  # Close all existing figures
    fig, ax = plt.subplots()
    markerline, stemlines, baseline = ax.stem(x_values, h_data, '-.')
    plt.setp(baseline, 'color', 'k', 'linewidth', 2)
    plt.title('Stem plot of selected filter')
    plt.xlabel('Index')
    plt.ylabel('Amplitude')

 

# Setup the widgets without any plotting
dropdown = widgets.Dropdown(options=list(h_variants.keys()), value='h32', description='Filter:')
output = widgets.Output()

# Function that updates the plot based on dropdown
def update_plot(change):
    with output:
        clear_output(wait=True)  # Clears the previous plot
        plot_stem(change['new'])  # Calls the plotting function

# Observe dropdown for changes
dropdown.observe(update_plot, names='value')

# Display the dropdown and output widgets in the desired order
display(dropdown, output)

# Call the update_plot function initially if needed
update_plot({'new': dropdown.value})

import numpy as np
import plotly.graph_objs as go
from ipywidgets import interact, Dropdown



h = np.real(np.fft.ifft(H))
middle = int(len(h)/2)
h = np.hstack((h[middle:], h[:middle]))

h_variants = {
    'h32': h[middle-16:middle+16],
    'h64': h[middle-32:middle+32],
    'h128': h[middle-64:middle+64],
    'h160': h[middle-80:middle+80],
    'h256': h[middle-128:middle+128],
}

# Function to create and update the stem plot based on selected filter
def create_update_stem_plot(h_data):
    x_values = np.arange(len(h_data))
    fig = go.Figure()
    
    # Add stems
    for x, y in zip(x_values, h_data):
        fig.add_trace(go.Scatter(x=[x, x], y=[0, y], mode='lines', line=dict(color='#1F77B4')))
    
    # Add markers
    fig.add_trace(go.Scatter(x=x_values, y=h_data, mode='markers', marker=dict(color='#1F77B4')))
    
    # Update layout
    fig.update_layout(
        title='Stem plot of selected filter',
        xaxis_title='Index',
        yaxis_title='Amplitude',
        template='plotly_white'
    )
    fig.show()

# Interactive widget setup
@interact
def plot_interactive(filter=Dropdown(options=h_variants.keys(), value='h32', description='Filter:')):
    h_data = h_variants[filter]
    create_update_stem_plot(h_data)

import plotly.graph_objs as go
from scipy import signal
import numpy as np

# Define filter sizes
filter_sizes = [32, 64, 128, 160, 256]

# Compute frequency responses
freq32, resp32 = signal.freqz(h32)
freq64, resp64 = signal.freqz(h64)
freq128, resp128 = signal.freqz(h128)
freq160, resp160 = signal.freqz(h160)
freq256, resp256 = signal.freqz(h256)

# Compute frequency responses
freqs = {}
resps = {}
for filt in [32, 64, 128, 160, 256]:
    freqs[f'h{filt}'], resps[f'h{filt}'] = signal.freqz(eval(f'h{filt}'))

# Create initial traces for the plot
traces = []
for filt in filter_sizes:
    traces.append(go.Scatter(
        x=0.5 * Fs * freqs[f'h{filt}'] / np.pi,
        y=np.abs(resps[f'h{filt}']),
        mode='lines',
        name=f'h{filt}',
        visible=False  # Initially set traces to be hidden
    ))

# Create layout for the plot
layout = go.Layout(
    title='Frequency Response of Selected Filters',
    xaxis_title='Frequency (Hz)',
    yaxis_title='Gain',
    yaxis_type='log',
    template='plotly_white',
    xaxis=dict(range=[0, 3000])  # Set x-axis range from 0 to 3000 Hz
)

# Define callback function to update visibility of traces
def update_visibility(selected_filters):
    visibility = [filt in selected_filters for filt in filter_sizes]
    return visibility

# Create dropdown menu for selecting filters
dropdown_menu = []
for filt in filter_sizes:
    dropdown_menu.append(
        dict(
            args=[{'visible': update_visibility([filt])}],
            label=f'h{filt}',
            method='update'
        )
    )

# Append 'All' option to dropdown menu
dropdown_menu.insert(0, dict(args=[{'visible': update_visibility(filter_sizes)}], label='All', method='update'))

# Update layout with dropdown menu
layout['updatemenus'] = [
    dict(
        buttons=dropdown_menu,
        direction='down',
        pad={'r': 10, 't': 10},
        showactive=True,
        x=0.1,
        xanchor='left',
        y=1.1,
        yanchor='top'
    )
]

# Create figure
fig = go.Figure(data=traces, layout=layout)

# Show plot
fig.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[18], line 9
      6 filter_sizes = [32, 64, 128, 160, 256]
      8 # Compute frequency responses
----> 9 freq32, resp32 = signal.freqz(h32)
     10 freq64, resp64 = signal.freqz(h64)
     11 freq128, resp128 = signal.freqz(h128)

NameError: name 'h32' is not defined

Παράθυρα Hamming και Kaiser

Show code cell source Hide code cell source

w_hamming=signal.hamming(len(h64))
h64_hamming = np.multiply(h64,w_hamming)
w_kaiser=signal.kaiser(len(h64),5)
h64_kaiser = np.multiply(h64,w_kaiser)

freq,resp64_hamming = signal.freqz(h64_hamming)
freq,resp64_kaiser = signal.freqz(h64_kaiser)

fig = go.Figure()

# Add trace for h64 (rectangular window)
fig.add_trace(go.Scatter(x=0.5 * Fs * freq / np.pi, y=np.abs(resp64), mode='lines', name='Rectangular', line=dict(color='#1F77B4')))

# Add trace for h64 with Hamming window
fig.add_trace(go.Scatter(x=0.5 * Fs * freq / np.pi, y=np.abs(resp64_hamming), mode='lines', name='Hamming', line=dict(color='green')))

# Add trace for h64 with Kaiser window
fig.add_trace(go.Scatter(x=0.5 * Fs * freq / np.pi, y=np.abs(resp64_kaiser), mode='lines', name='Kaiser', line=dict(color='red')))

# Update layout for a semilog plot
fig.update_layout(
    title='Απόκριση συχνότητας βαθυπερατού φίλτρου<br>παράθυρα hamming (πράσινο) και kaiser (κόκκινο)<br>(το αντίστοιχο ορθογωνικό σε μπλε)', title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Κέρδος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True)
)

# Show figure
fig.show(renderer='notebook')
# Convert figure to HTML
fig_html = pio.to_html(fig, full_html=False, include_plotlyjs='cdn')
display(HTML(fig_html))

Φίλτρα ισοϋψών κυματώσεων

Show code cell source Hide code cell source

from dash import Dash, dcc, html
from dash.dependencies import Input, Output
import plotly.graph_objs as go
import numpy as np
from scipy import signal


# Define the filters using the Remez algorithm
filters = {
    'Equiripple Filter 32+1': signal.remez(32, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'Equiripple Filter 64+1': signal.remez(64, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'Equiripple Filter 128+1': signal.remez(128, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'Equiripple Filter 160+1': signal.remez(160, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'Equiripple Filter 256+1': signal.remez(256, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
}

app23 = Dash(__name__)

app23.layout = html.Div([
    # Title
    html.H4('Select Equiripple Filter Length:', style={
        'textAlign': 'center',
        'marginBottom': '10px',
        'padding': '10px',  # Add some padding around the text
        'borderRadius': '5px',  # Optional: round the corners of the background
        'fontSize': '24px'  # Increase font size
    }),
    
    # Checklist for selecting filters
    dcc.Checklist(
        id='filter-checklist',
        options=[{'label': name, 'value': name} for name in filters.keys()],
        value=[],  # No default selection
        labelStyle={'display': 'block'},  # Display checkboxes in block style for better readability
        inputStyle={'display': 'block',"marginRight": "10px", "transform": "scale(1.5)"},  # Make checkboxes slightly larger
    ),
    
    # Graph for displaying frequency response
    dcc.Graph(id='frequency-response-plot'),
])

@app23.callback(
    Output('frequency-response-plot', 'figure'),
    [Input('filter-checklist', 'value')]
)
def update_frequency_response(selected_filters):
    fig = go.Figure()
    
    for filter_name in selected_filters:
        filter_coeffs = filters[filter_name]
        freq, response = signal.freqz(filter_coeffs, fs=Fs)
        fig.add_trace(go.Scatter(x=freq, y=20 * np.log10(np.abs(response)), mode='lines', name=filter_name))
    
    fig.update_layout(
        title='Frequency Response of Selected Equiripple Filters',
        xaxis_title='Frequency (Hz)',
        yaxis_title='Gain (dB)',
        template='plotly_white'
    )
    return fig

if __name__ == '__main__':
    app23.run_server(mode='inline',debug=True, port=8073)

Εφαρμογή του φίλτρου

Show code cell source Hide code cell source

from dash import Dash, dcc, html
from dash.dependencies import Input, Output
import plotly.graph_objs as go
import numpy as np
from scipy import signal
from scipy.io.wavfile import write


# Design the filters (replace with your actual filter design)
lp_filters = {
    'lpass32': signal.remez(32, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'lpass64': signal.remez(64, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'lpass128': signal.remez(128, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'lpass160': signal.remez(160, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs),
    'lpass256': signal.remez(256, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs)
} 

for name, lp_filter_coeffs in lp_filters.items():
    # Apply the filter to the signal
    s_lp_filtered = signal.convolve(s, lp_filter_coeffs, mode='same') / np.sum(lp_filter_coeffs)

    # Generate a filename based on the filter name
    filename = f'assets/{name}.wav'

    # Normalize the filtered signal to prevent clipping
    s_lp_filtered_normalized = np.int16((s_lp_filtered / np.max(np.abs(s_lp_filtered))) * 32767)

    # Save the filtered signal as a WAV file
    write(filename, Fs, s_lp_filtered_normalized)

app24 = Dash(__name__)

app24.layout = html.Div([
    dcc.Dropdown(
        id='filter-dropdown',
        options=[{'label': name, 'value': name} for name in lp_filters.keys()],
        value='lpass64'
    ),
    dcc.Graph(id='psd-plot'),
    html.Audio(id='audio-player', controls=True, src='assets/lpass32.wav')
])

@app24.callback(
    Output('psd-plot', 'figure'),
    [Input('filter-dropdown', 'value')]
)
def update_psd_plot(selected_lp_filter_name):
    selected_lp_filter = lp_filters[selected_lp_filter_name]
    s_lp_filtered = signal.convolve(s, selected_lp_filter, mode='same') / np.sum(selected_lp_filter)
    f, Pxx_den_lp = signal.welch(s_lp_filtered, Fs, noverlap=128, nperseg=256)
    
    fig = go.Figure()
    fig.add_trace(go.Scatter(x=f, y=Pxx_den_lp, mode='lines', line=dict(color='#1F77B4')))
    fig.update_layout(
        title=f'Φάσμα φιλτραρισμένου σήματος: {selected_lp_filter_name}',
        xaxis_title='Συχνότητα (Hz)',
        yaxis_title='Πυκνότητα φάσματος ισχύος',
        yaxis_type='log',
        template='plotly_white'
    )
    return fig

@app24.callback(
    Output('audio-player', 'src'),
    [Input('filter-dropdown', 'value')]
)
def update_audio_source(selected_lp_filter_name):
    return f'/assets/{selected_lp_filter_name}.wav'

if __name__ == '__main__':
    app24.run_server(mode='inline',debug=True, port=8074)

Show code cell source Hide code cell source

# Παράμετροι σήματος
Fs = 8192  # Συχνότητα δειγματοληψίας
t = np.arange(0, 1.0, 1/Fs)  # Διάρκεια σήματος 1 δευτερόλεπτο

# Δημιουργία του σήματος
s1 = np.sin(2*np.pi*700*t) + np.sin(2*np.pi*900*t) + np.sin(2*np.pi*1400*t) + np.sin(2*np.pi*2500*t)

# Σχεδιασμός των φίλτρων Parks-McClellan
lpass32 = signal.remez(32, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs)
lpass64 = signal.remez(64, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs)
lpass128 = signal.remez(128, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs)
lpass160 = signal.remez(160, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs)
lpass256 = signal.remez(256, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs)

# Φιλτράρισμα του σήματος
s1_filtered32 = signal.lfilter(lpass32, 1.0, s1)
s1_filtered64 = signal.lfilter(lpass64, 1.0, s1)
s1_filtered128 = signal.lfilter(lpass128, 1.0, s1)
s1_filtered160 = signal.lfilter(lpass160, 1.0, s1)
s1_filtered256 = signal.lfilter(lpass256, 1.0, s1)

# Φασματική πυκνότητα των φιλτραρισμένων σημάτων
f, Pxx1_den = signal.welch(s1, Fs, nperseg=1024)
f32, Pxx1_den32 = signal.welch(s1_filtered32, Fs, nperseg=1024)
f64, Pxx1_den64 = signal.welch(s1_filtered64, Fs, nperseg=1024)
f128, Pxx1_den128 = signal.welch(s1_filtered128, Fs, nperseg=1024)
f160, Pxx1_den160 = signal.welch(s1_filtered160, Fs, nperseg=1024)
f256, Pxx1_den256 = signal.welch(s1_filtered256, Fs, nperseg=1024)

# Σχεδίαση των φασματικών πυκνοτήτων σε Plotly
fig = go.Figure()


# Αρχικό σήμα
fig.add_trace(go.Scatter(x=f, y=10*np.log10(Pxx1_den), mode='lines', name='Original Signal'))

# Φιλτραρισμένο σήμα 32
fig.add_trace(go.Scatter(x=f32, y=10*np.log10(Pxx1_den32), mode='lines', name='Filtered Signal 32'))

# Φιλτραρισμένο σήμα 64
fig.add_trace(go.Scatter(x=f64, y=10*np.log10(Pxx1_den64), mode='lines', name='Filtered Signal 64'))

# Φιλτραρισμένο σήμα 128
fig.add_trace(go.Scatter(x=f128, y=10*np.log10(Pxx1_den128), mode='lines', name='Filtered Signal 128'))

# Φιλτραρισμένο σήμα 160
fig.add_trace(go.Scatter(x=f160, y=10*np.log10(Pxx1_den160), mode='lines', name='Filtered Signal 256'))

# Ενημέρωση layout
fig.update_layout(title='Φασματική Πυκνότητα Ισχύος',
                  xaxis_title='Συχνότητα (Hz)',
                  yaxis_title='Πυκνότητα Ισχύος (dB)',
                  template='plotly_white')

fig.show(renderer='notebook')
# Convert figure to HTML
fig_html = pio.to_html(fig, full_html=False, include_plotlyjs='cdn')
display(HTML(fig_html))

Ζωνοπερατά φίλτρα

Με αναλυτικό υπολογισμό της κρουστικής απόκρισης και παράθυρο

Show code cell source Hide code cell source

# Με αναλυτικό υπολογισμό της κρουστικής απόκρισης και παράθυρο kaiser
f1=700; f2=1500;  
Ts=1/Fs
f2m1=(f2-f1); f2p1=(f2+f1)/2; N=256
t=np.arange(-(N-1),N-1,2)*Ts/2
hbp=2/Fs*np.divide(np.multiply(np.cos(2*np.pi*f2p1*t),np.sin(np.pi*f2m1*t))/np.pi,t)
hbpw=np.multiply(hbp,signal.kaiser(len(hbp),5))

s1_bp=signal.convolve(s1,hbp,'same')

Show code cell source Hide code cell source

f0, Pxx1_den = signal.welch(s1_bp, Fs, noverlap=128, nperseg=256)

# Create Plotly figure
fig = go.Figure()

# Add trace for the power spectral density
fig.add_trace(go.Scatter(x=f0, y=Pxx1_den, mode='lines', line=dict(color='blue')))

# Update layout for a semilog plot
fig.update_layout(
    title='Φάσμα ζωνοπερατού σήματος sonar',title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Πυκνότητα φάσματος ισχύος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True, range=[np.log10(1e-15), np.log10(1e-7)])  # Adjust y-axis range
)

# Show figure
fig.show(renderer='notebook')
# Convert figure to HTML
fig_html = pio.to_html(fig, full_html=False, include_plotlyjs='cdn')
display(HTML(fig_html))

def play_sound(b):
    sd.play(20 * s1_bp, Fs)

# Create a button widget
play_button = widgets.Button(description="Play Sound")

# Link the button to the function that plays the sound
play_button.on_click(play_sound)

# Display the button
display(play_button)

Ζωνοπερατό ισουψών κυματώσεων

Show code cell source Hide code cell source

bpass = signal.remez(128, [0, f1*0.95, f1*1.1, f2*0.95, f2*1.05, Fs/2], [0, 1, 0], fs=Fs)
freq,resp_pm = signal.freqz(bpass)

# Compute frequency response for the equiripple bandpass filter
freq, resp_pm = signal.freqz(bpass)

# Create Plotly figure
fig = go.Figure()

# Add trace for equiripple bandpass filter response
fig.add_trace(go.Scatter(x=0.5 * Fs * freq / np.pi, y=np.abs(resp_pm), mode='lines', name='Equiripple Bandpass', line=dict(color='green')))

# Update layout for a semilog plot
fig.update_layout(
    title='Απόκριση συχνότητας ζωνοπερατού φίλτρου equirriple',title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Κέρδος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True)
)

# Show figure
fig.show(renderer='notebook')

# Convert figure to HTML
fig_html = pio.to_html(fig, full_html=False, include_plotlyjs='cdn')
display(HTML(fig_html))